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Abstract 

We study with Monte Carlo methods an ensemble of c = —5 gravity graphs, 
generated by coupling a conformal field theory with central charge c = —5 to 
two-dimensional quantum gravity. We measure the fractal properties of the 
ensemble, such as the string susceptibility exponent 7^ and the intrinsic fractal 
dimensions d/f. We find 7s — —1.5(1) and dn = 3.36(4), in reasonable agree- 
ment with theoretical predictions. In addition, we study the critical behavior of 
an Ising model on a quenched ensemble of the c = — 5 graphs and show that it 
agrees, within numerical accuracy, with theoretical predictions for the critical 
behavior of an Ising model coupled dynamically to two-dimensional quantum 
gravity, provided the total central charge of the matter sector is c = —5. From 
this we conjecture that the critical behavior of the Ising model is determined 
solely by the average fractal properties of the graphs, the coupling to the ge- 
ometry not playing an important role. 



1 Introduction 

Randomness in statistical systems arises in a variety of situations and is a very rich 
and complex subject. Quenched randomness is frequently used in studying the role 
of impurities and inhomogeneities in real physical systems where the characteris- 
tic time-scale of the disorder is much longer than other dynamics of the system. 
Annealed randomness, on the other hand, arises naturally in studies of fluctuating 
geometries, such as two-dimensional quantum gravity or fluid membranes, where the 
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disorder is dynamically modified by interaction between the geometry and matter 
fields living on the surfaces. 

For a statistical system coupled to random disorder, either in a quenched or 
annealed approach, the main question is to assess the effect randomness has on the 
critical behavior of the pure system. One prediction in this direction is the Harris 
conjecture |j^] which states that randomness changes the values of critical exponents 
only if the specific heat exponent a of the pure system is positive. This conjecture 
has been studied in many models with quenched disorder, such as the 2d Ising model 
|§] (where the Harris criterion is ambiguous as a = 0) and the Potts model For 
both models a change in the critical behavior is observed. 

All the above mentioned studies deal with weak disorder. More recently the 
critical behavior of systems on lattices with fractal structure very different from a flat 
surface has been investigated. Such systems arise naturally when matter, in the form 
of conformal field theories, is coupled to two-dimensional quantum gravity. These 
models can be studied either in a continuum formulation, by Liouville field theory, or 
using discretized approaches like, for example, models of dynamical triangulations, 
formulated either as matrix models or studied with numerical simulations. For these 
systems the disorder is, however, different from the one discussed above in that it is 
annealed, i.e. the models couple dynamically to fluctuations in the geometry. 

A remarkable degree of universality does emerge for models coupled to two- 
dimensional quantum gravity. Namely, the change in the critical behavior of the 
systems, and their back-reaction on the geometry, only depends on the total central 
charge of the matter sector. This manifests itself in the so-called KPZ scaling rela- 
tion which describe how the scaling dimensions of conformal operators are changed 
by the interaction with gravity Moreover, this universality also extents to the 
fractal structure of the surfaces, from which we derive the string susceptibility ex- 
ponent 7s and the fractal dimension dn- 

In view of this universality it is tempting to conjecture that the critical be- 
havior of a particular system, when coupled to a fluctuating geometry, only depends 
on the (average) fractal structure of the surface. Details of the interaction between 
the system and the geometry, or the geometrical fluctuations, are not important 
as such — they only serve the purpose of defining the average fractal geometry. If 
this conjecture is true it implies that how the average over disorder is performed, 
i.e. that the disorder is annealed, is not essential. In particular, predictions of the 
KPZ scaling relation for the change in the critical behavior should just as well apply 
to models with quenched disorder, provided the quenched average is taken over the 
same ensemble of disorder as is generated in the annealed approach. 

There are some recent simulations that have addressed the question of the 
critical behavior of spin models on a quenched ensemble of graphs generated by 
two-dimensional quantum gravity. Both the Ising model (H and the 10-state Potts 
model Q have been studied on an ensemble of pure gravity graphs (c = 0). For the 
Ising model a critical behavior compatible with an Ising model coupled dynamically 
to gravity was found, although the accuracy of the results is not sufficient to rule 
out the conjecture discussed above. 

The goal of this paper is two-fold. First, we want to investigate the fractal 
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geometry of two-dimensional quantum gravity coupled to a conformal field theory 
with central charge c = —5. More precisely, we want to determine the fractal di- 
mension of the corresponding surfaces, using recently developed finite-size scaling 
methods P, ^ and to compare it to the (contradictory) theoretical predictions that 
exist |lO|. Second, we want to investigate the critical behavior of an Ising model 
on a quenched ensemble of c = — 5 graphs and to compare it with predictions from 
Liouville theory, for the critical behavior of an Ising model coupled dynamically to 
two-dimensional quantum gravity, and to verify, or disprove, our conjecture about 
the effect of the disorder. Our motivation for choosing c = —5 is that both its pre- 
dicted fractal structure and the critical behavior of the Ising model is substantially 
different from both a flat space and for a pure two-dimensional quantum gravity. 
This makes these different critical behavior easier to distinguish in numerical simu- 
lations. 

The paper is organized as follows: In Section 2 we study the fractal properties 
of a c = —5 conformal field theory coupled to two-dimensional gravity. We define 
the model in Section 2.1 and discuss the details of the simulations in Section 2.2. 
In Sections 2.3 and 2.4 we present our measurements of the string susceptibility 
exponent 7^ and of the fractal dimension du- And in Section 2.5 we comment on 
how this particular ensemble of graphs differs from other types of graphs frequently 
used in studying disordered system. The second part of the paper deals with an 
Ising model on the c = —5 graphs in a quenched approach. In Section 3.1 we 
discuss the prediction from Liouville theory for the critical behavior of an Ising 
model coupled dynamically to two-dimensional gravity. In Section 3.2 we discuss 
details of the simulations and the observables we use to probe the critical behavior. 
In Sections 3.3 and 3.4 we determine the critical temperature of the Ising model 
and the corresponding critical exponents. Finally, in Section 4 we summarize and 
discuss our results. 



2 Geometrical properties of c = —5 graphs 

2.1 Dynamical triangulations coupled to scalar matter fields 

The model we study, and use to define our ensemble of random surfaces, is a dis- 
cretization of a free bosonic string theory embedded in D-dimensions (for an excel- 
lent review see e.g. Ref. |ll|). In the continuum formulation (Liouville theory) the 
partition function of two-dimensional quantum gravity coupled to a conformal field 
theory, with central charge c, is defined as 

Z{fl) = J P[5]P</.e-^/^'^^^-^A^('^'5) , (1) 

where p, is the cosmological constant, the integration is over equivalence classes of 
metrics [g], and 5'm(0, g) is the matter Lagrangian. For a free bosonic string, embed- 
ded in Z'-dimensions, Sm = / ^^iy/gg^^^daXd^jX , where X are the embedding 
coordinates of the string. 

Discretizing the model Eq. (|l|) using the simplicial gravity approach 1 12 ] , alias 



dynamical triangulations, the integration over metrics is replaced by a summation 
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over all triangulations T constructed by gluing equilateral triangles together along 
their edges into a simplicial manifold of a given topology [12|. The discretized 
(grand-canonical) partition function becomes 




where A is the area of the triangulation (number of vertices), x is the embedding of 
a vertex in a D-dimensional space, and (ij) indicates that the sum is over adjacent 
vertices in the triangulation. The center of mass is kept fixed to eliminate the 
translational zero mode. 

The integration over the Gaussian fields in Eq. (^) can be carried out explicitly 
and the canonical (fixed area) partition function becomes 

Za = J2 (d^t Cr)-^/2 , (3) 
rer(A) 

where Ct is the adjacency matrix of the the triangulation T: 

{Qi if « = j, 

—Cij if i and j are adjacent, (4) 

otherwise. 

Here qi is the order of vertex i and Cij is the number of edges connecting the adjacent 
vertices i and j. Note that as we use degenerate triangulations, which are defined 
below, there can be more than two edges connecting the vertices i and j and the 
order of vertices is defined excluding self-loops. In calculating the determinant of 
Ct one vertex is excluded from the graph in order to eliminate the zero mode. 

In the partition function Eq. (|3|) the embedding dimension D now appears as 
a free parameter and is no longer restricted to positive integer values. In particular, 
it can be taken to be negative; this corresponds to coupling a conformal field theory 
with a negative central charge c = D to two-dimensional gravity. 

The sum over triangulations T in Eq. (|2|) is over an appropriate class of trian- 
gulations T of area A and with fixed topology. Different classes amount to different 
discretization of the manifolds. By universality arguments different choices of T 
should yield the same critical behavior as long as they only differ at the level of the 
discretization. This statement is known to be true for two commonly used classes 
of triangulations, degenerate (Tb) and combinatorial (Tc) triangulations. Combi- 
natorial triangulations are defined by forbidding self-loops (an edge starting and 
ending at the same vertex) and vertices connected by more than one edge. These 
pathologies are, on the other hand, allowed in the class of degenerate triangulations. 
For these particular classes of triangulations the model Eq. (^ has been solved for 
some special values of c This universality holds even if curvature fluctuations 
of the triangulation, i.e. the vertex orders qi, are maximally restricted, to = 6ibl, 
as shown by numerical simulations 

Finite-size effects do, however, depend on the discretization; in particular, for 
degenerate triangulations they are substantially smaller (one order of magnitude) 
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Table 1: The number Ng of c 



-5 graphs, of area A, we have generated by simulating the 
model Eq. (^. Also shown is the number Nr of different replicas (graphs), for each area, 
on which we simulated the Ising model (Section 3). 



A 

100 
150 
200 
300 
400 
600 
800 
1200 
1600 



3000 
3000 
3000 
3000 
2000 
2500 
700 
196 
96 



Nr 

200 

200 

130 

60 
178 
96 



than for combinatorial |15|. In the work presented in this paper we use degenerate 
triangulations of spherical topology. 

It is possible to get some insight into the behavior of the model in the limits 
D ^ iboo by looking at what triangulations dominate the sum in Eq. (^) As the 
determinant of Ct is related to the number of spanning trees on a particular graph, 
the graphs dominating in the limit D +oo are those with a minimal number 
of spanning trees. This corresponds to triangulations with a branched polymer 
structure. In the limit c — > — oo, on the other hand, we expect the determinant Ct 
to take its largest value for triangulations with a maximal number of spanning trees. 
This implies that the dominant triangulations in this limit will be flat with qi = 6 
(apart from few defects Qi ^ 6 needed to form a closed spherical surface). 

2.2 Simulating a c = —5 gravity theory 

We have studied numerically the partition function Eq. (^) for c = — 5 using Monte 
Carlo simulations. The space of all triangulations is explored using the so-called 
edge- flip algorithm in which an edge kj, common to two triangles tijh and tkji, 
is removed (or flipped) and replaced by the edge //j^. This algorithm is known to be 
ergodic (for fixed area). Unfortunately, the determinant in Eq. (^) corresponds to a 
non-local action making the simulations very difficult. Every time a flip of an edge 
is proposed the whole determinant has to be recalculated, an operation requiring 
on the order of N"^ floating point operations. We can exploit the locality of the flip 
and using standard techniques to reduce this to A^^ floating point operations 
per flip, but this still limited our simulations to triangulations consisting of up to 
1600 vertices. Nevertheless, those triangulations are considerably larger than have 



been simulated before with this method [18|. 



In the simulations we stored triangulations after every 10 to 20 Monte Carlo 
sweeps, were each sweep consisted of flipping, approximately, an area worth of edges. 
With the auto-correlations present, which are relatively modest on such small sur- 
faces, the stored graphs were more or less independent. In Table 1 we show the total 
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number of graphs we generated, and analyzed, for each area. 



2.3 The string susceptibility exponent 7^ 

To determine the critical behavior of the model Eq. (^), with c = —5, we mea- 
sured two critical exponents: the string susceptibility exponent 7^ and the fractal or 
Hausdorff dimension d^. The latter will be discussed in the next section, here we 
consider 7<j. 

The string susceptibility exponent is defined by the singular behavior of the 
grand-canonical partition function Eq. (Q) as the cosmological constant fi approaches 
its critical value: 

Z{fl) ^ Zreg + {fi-ficf-^\ (5) 

This implies that the canonical partition function behaves asymptotically as 

Z{A) ~ e l^cAj^-r.-s ^ ^ ^ ^ (g) 

The value of 7s can be calculated from the Liouville field theory Q: 

7s = ^ (c- 1 - V(c-25)(c-l)) . (7) 

For c > 1 this implies a complex critical exponent and the corresponding theory is 
not well defined — this is related to the existence of tachyons in string theories in 
embedding dimension larger than two. For c < 1, on the other hand, the theory is 
well defined and, specifically, for c = — 5 Eq. (0) predicts 7^ = —1.618 . . . 

In order to measure 7<j for the ensemble of triangulations, generated in our 
simulations, we study the distribution of so-called minhus (minimal neck baby uni- 



verses) [19|. A minbu is a part of the triangulation connected to the rest through a 
minimal neck; for a degenerate triangulation a minimal neck is simply a self-loop. 
By counting in how many ways a minbu of size B can be connected to a surface of 
size {A — B), the size distribution of minbus can be written as: 

~ [{A - B) ^ (g) 

In the last step we have used the asymptotic behavior of the partition function, 
Eq. (|6|) — in this way the leading exponential behavior cancels out and we can 
directly access the sub- leading corrections governed by 7^. 

The distribution n^iB) is easily measured in numerical simulations and 7^ 
extracted by a fit to Eq. (P). In practice, though, one has to impose a lower cut-off 
on the size of minbus included in the fit as there are finite-size corrections to the 
asymptotic form Eq. (P). This can be done in a systematic way; small minbus are 
thrown away in the fitting procedure until one gets a stable value of 7^ and an 



acceptable x value for the fit [20|. Unfortunately, the large negative value of 7^ 
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for the c = — 5 graphs makes this measurements more difficult than for previously 
studied ensembles of triangulations, as the distribution falls off very rapidly. 

We have measured 7^ for surfaces of area up to 300 vertices; for the larger 
surfaces our statistics was not sufficient for a reliable determination of the minbu 
distribution. Although these are very small surfaces, our experience in measuring 
7s for other models of dynamical triangulations shows that they are large enough, 



provided one is using degenerate triangulations |15|. In this particular case the lack 
of statistics is much more of a problem than the smallness of the triangulations. But 
as determining 7^ is of secondary interest for us, we did not deem it worthwhile to 
invest too much CPU-power in increasing the statistics. We get the following values 
of 7s : 



A 


7s 


100 


-1.52(17) 


150 


-1.42(15) 


200 


-1.37(22) 


300 


-1.71(19) 



Those values are in reasonable agreement with the theoretical value 7^ 



-1.61^ 



2.4 The vertex-vertex correlation function and dn 

Another exponent that characterizes the fractal structure of the surfaces is the 
(intrinsic) fractal dimension dn- It can be defined from the behavior of the vertex- 
vertex (or two-point) correlation function QAi''')'- 

9A{r) = j{^6{d,,-r))T, (10) 

where dij is the (minimal) geodesic distance between two vertices i and j and the 
statistical average is performed over all triangulations T. This correlation function 
simply counts the number of vertices (or the area of the manifold) at a geodesic 
distance r from a marked vertex i, averaged over all vertices i. We expect its short 
distance behavior to be QAif) ~ r'^^"^, provided r <C A^^'^^ . 

On a triangulation we define the geodesic distance between two vertices as 
the shortest path between them traversed along links. Alternatively, one can define 
a triangle-triangle correlation function tA{r), analogous to Eq. (|lO|), in which case 
the geodesic distance is defined as the shortest path between two triangles traversed 
through the center of triangles. Although those two definitions will result in very 
different correlation functions for a particular triangulation, they should define the 
same fractal dimension. 

To extract the fractal dimension from the measurements of the correlation 
function Eq. ( [lO| ) it is convenient to use methods of finite-size scaling. Assuming 
that the only relevant length-scale in the model is defined by A^^'^^ , general scaling 
arguments ^ |8| imply that 

9A{r) ~ A'-'/^^ Fix), (11) 
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where we have introduced the scahng variable 



x=—r7i-- (12) 

In fixing the prefactor in Eq. (|l^) we have used that ^rdAij) = A. 

In practice there are strong finite-size corrections to this scaling behavior. 
The measurements of dn can be improved, considerably, by including the simplest 
finite-size correction to x, i.e. by introducing the so-called shift j^: 

where the shift parameter a will in general depend on the particular observable we 
consider. This scaling correction has been applied successfully to all geometric and 
matter correlation functions, even with quite different functional dependence on x 



|21|. Its possible geometric origin has been investigated analytically for c = in 
|22|. The shift was also found and calculated analytically in correlation functions on 
branched polymers where it has been shown to contain the singularity responsible 
for the crumpling phase transition ||23|| . 

It is possible to generalize Eq. ( p!o| ) to higher moments of the correlation func- 
tion. Introducing the loop- length distribution pyi(r, /), which counts the number of 
loops of length / at geodesic distance r from a marked vertex z, we define the A;-th 
moment of / as [10, 22| 

l\{r) = Y^l'^ pA{r,l) . (14) 
I 

Note that 5yi(r) = l\{r). On a triangulation loops are defined as a connected 
(by links) subset of the vertices that are at distance r from However, the 
generalization of the scaling hypothesis Eq. (|Tl|) is more subtle, as it depends 



on how the boundary length I scales with the radius r [22|. If we assume that 
dim[r^] = dim[i] = dim[^^/'^^], then we have for the moments: 

Z^(r) = ^2fc/rfii ^ ^ > i_ (15) 

Alternatively, if we assume that dim[/^] = dim[^] = dim[r'^^], we get: 

I'Xir) = A'^/^ Fu{x). (16) 



Our measurements clearly favor the former scaling form Eq. (15); this is compatible 
with results recently obtained for a c = —2 theory |^2|. This implies, from dimen- 
sional analysis, that the boundary length / acquires an anomalous scaling dimension 
with the area; at present the origin of this behavior is still not understood. Note 
that in the only soluble case, c = 0, both scaling forms agree. 

^In a similar way we can define a loop-length distribution for the triangles at a distance r from 
a marked triangle. In that case we define a loop as a connected subset of triangles sharing at least 
one vertex. If we define loops as triangles sharing at least one link, then those loops always stay at 
the level of the lattice cut-off. 
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Table 2: The fractal dimension dfj for an ensemble of c = —5 graphs, determined by 
collapsing the correlation functions onto a single curve using the scaling forms Eqs. ( pi] ) 
and and the shift Eq. (^3|). The corresponding (optimal) value of the shift parameter 
a is included. Graphs of area 200 to 1600 were included in the analysis. The corresponding 
values for the triangle-triangle correlation function tA{r) are also included. 



/^(r) k 




9A{r) 1 
2 
3 
4 

tA{r) 


3.36(16) 0.6(6) 
3.33(5) 0.4(3) 
3.36(4) 0.4(3) 
3.36(4) 0.4(3) 

3.07(24) 2.4(1.6) 



We have measured the moments l\{r) for k < A, together with the triangle- 



triangle correlation function tA_{r)- The scaling behavior, Eqs. (11) and (|T5|), is 
analyzed by "collapsing" the distributions corresponding to different area onto a 
single curve. For this purpose it is convenient to interpolate between the (discrete) 
values of X — this we do using a cubic-spline (see Ref . for details) . The collapse 
depends on only two free parameters, a and dn', their optimal values are determined 
by minimizing the x^^value for the collapse, where the x^~value is defined by the 
difference between the curves after rescaling. The results are shown in Table 2. We 
used triangulations of area 200 to 1600 in the analysis; if we included smaller graphs 
it was not possible to get an acceptable x^-value for the fit. We observe, as has 
been observed before |^, ^, that the finite-size effects are considerable bigger for the 
triangle-triangle correlation function. This is also reflected in much larger value of 
the shift parameter a. 

To demonstrate how good this scaling behavior actually is, we show in Figure 1 
the scaled curves for the vertex- vertex correlation function 17^1 (r). These curves have 
been scaled using dn = 3.36(4), the average value of the fractal dimension from 
Table 2, and the corresponding value of the shift. The quality of the scaling for the 
higher moments, and for tAir), is equally impressive. 

Finally, we look at the loop-length distribution pAif, 0- Although, in principle, 
it contains the complete set of informations about the moments, it is not very 
convenient for extracting the fractal dimension. Nevertheless, we demonstrate that 
its scaling behavior is consistent with the measured fractal dimension. Since our 
measurements show that for the higher moments / ~ (see Eq. (|l5[)), i.e. the same 
behavior as for c = 0, we expect that 

PA{r,l) « r^F,{l/r^). (17) 

We show an example of this scaling in Figure 2. The relation Eq. ([T7|) is, of course, 



exact only for Poo{r,l), as is the case for c = |1C], and finite-size effects are ex- 
pected for small Z/r^. The deviation from the above scaling behavior is, in principle, 
determined by the scaling corrections Eqs. (pT|) and (^) |2l|] . 
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Figure 1: The scaled vertex-vertex correlation function g'j^ir) = gA{r)/A^^^^'^'' vs. the 
scaling variable x = {r + a)/A^^'^" . The measured value of the fractal dimension dn = 3.36 
and a shift a — Q.Q were used in the scaling. 



As we mentioned in the Introduction, there exist several theoretical predictions 
for the fractal dimension for theories of two-dimensional quantum gravity coupled 
to matter. Using string field theory, or a transfer matrix approach, with a modified 
definition of a geodesic distance, one gets |^: 

2 24 

dn = ^ = = . (18) 

\ls\ 1-C+ v^(25-c)(l-c) 

An alternative expression is obtained by studying the diffusion on the surface, within 
the framework of Liouville theory [|lO| : 



, ^ x/25^ + V49^ 
These two analytical predictions disagree, except in the the case of pure gravity; 



the only case where an exact solution exist [24|. For c = we have du = 4, this 
value is also obtained in numerical simulations using the scaling methods described 
above 0, On the other hand, both Eqs. (18) and ( [l9|) disagree with the results of 



extensive numerical simulations of conformal field theories with < c < 1 coupled to 
gravity [0, ^, |2^, ^] . The numerical simulations indicate that the fractal dimension 
is 4 independent of the coupling to matter. In [^], the discrepancy for c = 1/2 
— one Ising model coupled to gravity — has been attributed to the subtlety in 
how the continuum limit should be taken is this model. An alternative proposal 
has been made in [^] where the authors suggest that a different scaling variable, 
related to the size of clusters, should be used in the derivation of Eq. (|l^) instead of 
the geodesic distance. It should be emphasized, however, that the present numerical 
accuracy is not enough to rule out completely the prediction of Eq. ([T9|) for c > 0. 
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Figure 2: An example of the scaling behavior of the loop-length distribution function 
PAir, I) for A = 1600. The distributions for different choices of I are re-scaled in accordance 
with Eq. 0. 



The only model with c < where the fractal dimension has been measured 
numerically is c = —2 or topological gravity [22|. This is a very special case were it 
is possible to sample the space of triangulations recursively; this allow simulations 
of very large graphs (graphs of up to 8 x 10^ triangles were used in Ref. |^^). The 
measured fractal dimension in this case is dn = 3.58(4), which agrees very well with 
the prediction from Liouville theory, Eq. (pj|): dn = 3.561 . . . 



Our result for c = —5, dn = 3.36(4), also agrees reasonably well with Eq. (|19|), 
which predicts dn = 3.236 ... for c = — 5, especially given the smallness of our 
graphs. On the other hand, it completely rules out the prediction of Eq. (|l8|): 
dn ~ 1.236. Although the value 3.36(4) is slightly larger than predicted by Eq. (p!9[), 
we notice that the fractal dimension obtained using the triangle-triangle correlation 
function: dn = 3.07(24), is somewhat smaller . This effect was also noticed for the 
c = — 2 model where the two values converged to Eq. (|l8|) on larger graphs. This is 
what we observe in our simulations as well. 



2.5 Comparison with other random lattices 

The fractal structure of the c = — 5 graphs differ substantially from other types 
of random lattices frequently studied. An ensemble of random lattices, commonly 
used in the study of quenched disorder, is Poissonian random lattices. They are con- 
structed by distributing vertices uniformly on a two-dimensional manifold and link 
them together to form a triangulation, usually following a prescription by Dirichlet 
and Voronoi [ p8[ . 

We have compared the properties of our ensemble of graphs to that of Voronoi 
triangulations by looking at the probability distribution of vertex orders pn- In 
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Figure 3: The (normalized) curvature distribution p„ for an ensemble of graphs correspond- 
ing to 2(i-gravity coupled to matter with c = — 5, and 5, respectively, and for Voronoi 
triangulations. The interpolating curves are just to guide the eye. 

Figure 3 we plot this distribution for the c = —5 graphs, Voronoi triangulations, 
and, for comparison, pure gravity (c = 0) and branched polymer (c = 5) graphs. 
As in all cases the distributions are for triangulations, they all have the mean value 
Pn = 6, but in other aspects they differ. The distribution for Voronoi graphs is 
peaked sharply around the mean value and falls off rapidly as n increases. The pure 
gravity distribution peaks at the smallest possible curvature (n = 1 for degenerate 
triangulations) and falls off much slower. The distribution for c = —5 lies in-between, 
as could be expected from its fractal dimension^ 

As Figure 3 indicates, the randomness of the Voronoi triangulations is only 
at the local level, their global fractal structure, defined by, for example, the fractal 
dimension or the string susceptibility exponent, will be the same as of the under- 
lying two-dimensional flat lattice used in constructing them. In this respect their 
randomness is only a (small) local perturbation around the flat background. This is 
also reflected in that the critical behavior of spin models on such triangulations is 
the same as for a flat lattice |2^. This is contrary to the other ensembles of graphs 
shown in Figure 3, which have a genuinely different global fractal structure. 

It should also be emphasized that although local observables such as pn reveal 
something about the randomness of the graphs, they are not a universal property and 
not a good indicator for the critical behavior of the corresponding model. Within 
the context of 2d gravity the local properties of a given model can be changed by 
adding an irrelevant operator to the action in Eq. (|^; for example, a term that 
couples to the local curvature. But, except in extreme cases, this modification does 

^It is worth noting that the fractal dimension alone is not enough to characterize the graphs. 
For example, both a flat lattice and branched polymers have dH ~ 2 
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not affect the critical behavior of the model |14, 30 1. 

It is, however, possible to change the global fractal properties of the graphs 
defined by the 2(i-gravity model Eq. (^) by adjusting the embedding dimension D. 
In particular, by taking D negative we can create an ensemble of triangulations with 
any fractal dimension between 2 and 4. This makes the model of Eq. (^) ideal for 
investigating how the fractal properties, like the fractal dimension, affect the critical 
behavior of a spin model living on the graphs. As will be demonstrated in the next 
section, in the case of the Ising model this results in critical exponents that are 
radically different from those of the Ising model on a flat lattice. 

3 The Ising model on c = — 5 graphs 
3.1 Predictions from the continuum 

In this second part of the paper we investigate how the critical behavior of an Ising 
model is modified when it is defined on a quenched ensemble of c = —5 graphs. In the 
annealed case coupling a conformal field theory to 2(i-gravity results in dressing the 
dimensions of primary operators of the theory. The dressed weights are given by the 
KPZ formula Q. From those weights one can calculate the new critical exponents. 
As a dynamical triangulation, the only case that has been solved explicitly is the one 



Ising model coupled to gravity, formulated as a two-matrix model |13|. In that case 
the phase transition changes from second to third order and the critical exponents 
agree with the ones calculated using the KPZ formula. 

All these calculations emphasized above, for the Ising model on an 

annealed ensemble of graphs. Now the conjecture, discussed in the Introduction, is 
that the interaction between the Ising model and the geometry is not important; 
the critical behavior is simply determined by the average fractal structure of the 
triangulations. This implies that we can also use the KPZ formula to calculate the 
critical behavior for an Ising model on a quenched ensemble of graphs, provided 
they define the appropriate fractal structure. If the conjecture is true the critical 
behavior of an Ising model on a quenched ensemble of c = —5 graphs is the same as 
that of the Ising model coupled dynamically to gravity with additional c = —11/2 
conformal field coupled to it. 

Given a primary operator of conformal weight A*^ in the original theory, the 
KPZ formula gives its conformal weight after coupling to gravity: 

For the Ising model the relevant operators are the energy density e and the spin 
cr; from the conformal weight of those operators we can calculate the specific heat 
exponent a and the magnetization exponent f3: 
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In flat space the conformal weights of the Ising operators are: = 1 and = |; 
hence the KPZ formula, together with Eq. (^T|), predicts a ~ —0.452432 and /3 ~ 
0.234186 for an Ising model on c = —5 graphs. The full list of exponents is shown in 
Table 5. As those exponents are very different from both the Ising model exponents 
on a flat lattice (the Onsager exponents) and coupled dynamically to pure gravity 
(both sets of exponents are included in the table) it should be relatively easy to 
distinguish them in numerical simulations. 



3.2 The Ising model simulations 

We have simulated the Ising model for several independent replicas of the c = —5 
graphs we generated. We did this for graphs of size 100, 200, 400, 800, 1200 and 
1600 vertices — for the number of replicas we sampled see Table 1. We use the 
standard definition of the Ising model: 

^ising(T) = ^e^^fe)^*^^' , (22) 

where Oi = ±1 and {ij) denotes adjacent vertices in the graph T. The Ising spins are 
placed on the vertices of the graphs; alternatively they could be put on the triangles, 
but based on duality arguments we would expect the same critical behavior (this is 
the case for an Ising model coupled dynamically to 2(i"gravity [|^]). Naively, one 
might expect that placing the spins on the triangles was preferable as for a graph 
of a given area this yield twice as many spins. But as the finite-size effects are 
dominated by the geometry not the number of spins this not the case. In fact, it 
would only increase the computational efforts in updating the spin configuration. 

We simulated at several values of the inverse temperature /? = J/ksT for 
P £ [0.15,0.25]. A Swendsen-Wang cluster algorithm was used to update the spin 
configuration and typically about 200.000 measurements taken at each /3 value; the 
measurements were separated by 10 to 20 cluster updates. For each measurement 
we stored the total energy of the system, E = Yl{ij) ^i^ji the magnetization, 
M = Yli^i'j other observables can be constructed from the those two. To in- 
terpolate between measurements at different temperatures we used standard multi- 



histogram methods |31|. 



All observables O were calculated for each replica independently and the av- 
erage over the different replicas r performed afterwards: 



^E^- (23) 



Nr 



This is in accordance with the philosophy that the quenched average should be 
performed at the level of the free energy, not at the level of the partition function 
|l32|| . The error on O is estimated by a jackknife analysis over different replicas. 

To analyze the critical behavior of the Ising model we constructed several 
standard observables. In addition to the energy density (e) and magnetization (m) 
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per vertex we considered: 

Ca = /3M((e2)-(e)2) 
A ((m^) - (m) >2) 

(e)4 



X 
Ua 



cq + ciA"'/'"^" (specific heat), 
j^-r/udH (magnetic suscept.), 



-0 



Va 



(Binder's cumulant), 
(energy cumulant), 



(24) 



together with the various derivatives of the magnetization: 



^\rn\ 
I-^ln \m\ 



d(|m|) 
d/3 

dln(|m|) 
dln(m^) 
AUa 



A{{e){m) - (e|m|)) 
A({e) 



(e|m|) 



^(1 - f/A) ((e>-2i^ + i 



(m4> 



^(l-/3)/z/d^f 



(25) 



In Eqs. (|24| ) and ( P5[) we have included the expected finite-size scaling behavior. 
For some of the observables this scaling behavior applies both at the infinite-area 
critical temperature (3c and to the scaling of their extremal values. 

3.3 The critical temperature 

To determine the critical temperature (3c we use that the pseudo-critical temperature 
(3c{A), defined by the location of peaks in the different observables, approaches the 
infinite-area value like: 



MA) - Pc + 



Co 



(26) 



The observables we have used for this purpose are: Ca, ^\m,\ > |m| > ^\nm? > 
and D;7^, all of which have have well resolved peaks. An example of this is shown in 
Figure 4 were we plot the specific heat and the magnetic susceptibility. In the figure 
the curves shown are the averages over different replicas of the model, Eq. (23). 
To determine (3c{A)^ and to estimate the corresponding error, we locate the peaks 
on each replica independently and then take the averages over the replicas. In 
the infinite-area limit, were one expects the triangulations to become self-averaging 
with respect to the fractal structure, the distribution of (31{A) should approach a 
delta- function centered at (3c- This we demonstrate in Figure 5 where we show the 
(normalized) distributions of pseudo-critical temperatures p{(3c)', in this case for the 
specific heat. And, indeed, the distributions get narrower as the area is increased. 
Similar behavior is observed for the other observables. 

The fit to Eq. (^) is made considerably easier by an independent determina- 
tion of the critical exponent udn- The peak values of the derivatives Dinimh ^inm^, 
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I 

0.16 0.18 0.2 0.22 0.24 



Figure 4: (a) The specific licat Ca for an Ising model on an ensemble of quenched c = —5 
graphs. (&) The corresponding magnetic susceptibility For both observables the curves, 
for each area A, are the averages over different replicas and the errors (dashed lines) are 
estimated using jack-knifing over replicas. 
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Figure 5: Fits of the (normalized) distributions of the pseudo-critical temperature /3'^{A) 
to a Gaussian distribution: p{/3^) — a exp(— (/3^ — /b). This is for /3c (^) defined by 
peaks in the specific heat C^, for an Ising model simulated on different replicas r of c = —5 
graphs. Note that the symbols are only to identify distributions corresponding to different 
area. 




0.04 0.08 0.12 0.16 -1/,^ 

Figure 6: Scaling of the pseudo-critical temperature i3c{A) vs. the scaled area A^^/'^d-n ^ 
using vdH — 2.464. This is shown for all the different observables we considered. The 
approach to the (infinite-area) critical temperature /3c is demonstrated by a linear fit to 
Eq. ( p6| ) (excluding the smallest system size). 
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Table 3: The critical temperature (3c, for the Ising model on c — —5 graphs, determined 
from the scaling of the location of peaks in the different observables. This is both for 
Am = 100 and 200 as the smallest area include in the linear fit to Eq. (p6|). Also included 
is the x^-value (per d.o.f.) for the fits. 



o 


A e {100,1600} 

/3c 


A e {200, 1600} 

X^ 


Ca 


0.2153(15) 


19 


0.2140(17) 7.9 


X 


0.2197(11) 


6.7 


0.2197(16) 7.9 


Ua 


0.2233(12) 


5.6 


0.2214(10) 2.1 




0.2187(5) 


1.3 


0.2183(7) 0.8 


^In m'^ 


0.2188(7) 


4.5 


0.2186(10) 4.0 




0.2184(5) 


0.3 


0.2185(4) 0.1 




0.2188(22) 


15 


0.2200(13) 2.0 


Average 


0.2188(9) 


0.2187(9) 



and D[/^ all scale like A^^'^'^" . We have used this to determine i/dn by a linear fit; 
those exponents are shown in Table ^ and the corresponding fits in Figure ^. To 
demonstrate the finite-size effects we have done the fits with both = 100 and 
Am = 200 as the smallest area included in the fit. Although the exponents are 
not substantially altered by excluding graphs of area 100, for most observables the 
X^-value of the fit is, however, unacceptably large if they are included. 

Using the average value udn = 2.464(19), we obtain the critical temperature 
from a linear fit to Eq. (^). The result, for the different observables we considered, 
is shown in Table ^. This yields the average value /3c = 0.2187(9) (using Am = 200). 

3.4 The critical exponents 

We then proceed to determine the critical exponents of the model by fitting the 
different observables to their expected finite-size scaling behavior, Eqs. (|2^ and 
(p5|), both at the critical temperature and, were appropriate, also their peak values. 
The exponents, obtained from those fits, are collected in Table ^. 

The exponents have been determined both for Am = 200 and 400, except for 
the specific heat which is the most difficult quantity to analyze. It requires a non- 
linear 3-parameter fit to the scaling behavior Eq. ( p4|) and, in addition, as it does 
not diverge for this particular model the regular background term is all the more 
important. Thus we could not get a stable fit if we excluded graphs smaller that 
400 and, as is apparent from the quoted errors, the estimate of this exponent is the 
least reliable. On the other hand, for the other exponents, most determined from 
more that one observable, we get very consistent estimates. In all cases graphs of 
area 100 had to be excluded in order to obtain an acceptable fit. 

Although the triangulations we have employed in this study are smaller than 
those usually used in simulations of two-dimensional gravity — due to the difficulty 
in simulating the non-local action Eq. (^) — this is to a large extent compensated 
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Table 4: The critical exponents of the Ising model on a quenched ensemble of c = —5 
graphs, determined from the scaling of various observables both at the critical temperature 
~ 0.2187 and, where appropriate, of peak values. Graphs of area A > Am, with Am = 200 
and 400, are included in the fits to Eqs. (24) and ( p5|) (except for Ca where it was not possible 
to obtained a reliable fit if area 200 was excluded). 



o 


(A) Scaling of peaks 


(B) Scaling at f3c 


Am = 200 


Am = 400 


Am = 200 


Am = 400 


vdn An|m| 

I?l 2 

P/ydn \m\ 

D\m\ 

lIvdH X 
a/vdu Ca 


2.414(20) 7.8 
2.405(23) 9.4 
2.591(43) 2.4 

0.1001(69) 1.3 
0.7966(31) 2.9 
-0.266(90) 0.7 


2.469(37) 2.3 
2.456(23) 1.3 
2.528(72) 1.7 

0.1070(58) 0.1 
0.7922(99) 2.6 


2.413(15) 1.0 
2.399(22) 2.4 
2.50(12) 4.3 

0.1055(34) 3.0 
0.0967(52) 2.5 

0.7825(57) 6.8 

-0.269(78) 0.4 


2.451(28) 0.4 
2.431(35) 0.9 
2.45(11) 3.3 

0.1006(67) 0.5 
0.1064(79) 1.0 

0.7958(146) 4.0 



by a smaller fractal dimension, which sets the relevant length-scale that controls 
the finite-size effects. Hence we are able to obtain reliable estimates of the critical 
exponents, even from graphs of such modest size. 

We have collected in Table ^ our numerical estimates of the exponents of an 
Ising model on a quenched ensemble of c = —5 graphs. The exponents shown are 
a weighted average over the values in Table ^ (corresponding to Am = 400). For 
comparison, we also included in the table the critical exponents for the Ising model 
on a flat lattice (the Onsager exponents), for the Ising model coupled dynamically 
to gravity (c = 1/2), and for the Ising model coupled dynamically to gravity and 
conformal matter with central charge cm = —11/2 (or c = —5). Comparing those 
exponents it is clear that our result agrees very well with the last set of exponents. 
That the critical behavior agrees so well with the predictions from the KPZ scaling 
relation for c = — 5 matter coupled to gravity, strongly supports our conjecture about 
the effect of disorder on the critical behavior. That is, the dynamical interaction 
between the matter and the geometry is not important as such, only that they result 
in a well defined (average) fractal structure for the surfaces. 

4 Discussion 

The main results of the work presented in this paper can be summarized as follows: 

(a) The fractal dimension of surfaces, defined by a conformal field theory with 
central charge c = — 5 coupled to two-dimensional quantum gravity, is cIh = 
3.36(4). This is in reasonable agreement with, and supports, the theoretical 
prediction Eq. (|l9|), whereas it definitely rules out Eq. ([l8|) . 
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Table 5: The critical exponents of the Ising model on a flat two-dimensional lattice (the 
Onsager exponents), and coupled dynamically to 2(i-gravity both for c = 1/2 and c — —5. 
This is compared to the results of our simulations for the Ising model on a quenched ensemble 
of c = —5 graphs. Also shown are the corresponding exponents for the fractal structure of 
the graphs: 7s and dn- 





vdH 


a/vdu 


P/vdn 


l/vdn 


7s 


dn 


Onsager 


2 


log 


0.0625 


0.875 


— oo 


2 


c= -5 


2.452 ... 


-0.184 ... 


0.0955 ... 


0.809 ... 


-1.618 ... 


3.236 ... 


c= 1/2 


3 


-0.333 ... 


0.1666 ... 


0.666 ... 


-0.333 . . . 


4.0-4.2 


Quenched 


2.464(19) 


-0.27(8) 


0.105(4) 


0.793(8) 


-1.5(1) 


3.36(4) 



(6) The critical behavior of an Ising model on a quenched ensemble of c = —5 
graphs agrees well with the predictions, from the KPZ scaling relation, for an 
Ising model on an annealed ensemble of graphs with identical fractal properties. 

The first result, especially combined with the recent simulations of 2(i-gravity 
for c = —2 [22 1, lends a strong support to Eq. (|l^) as a correct description of 



the fractal structure of two-dimensional quantum gravity for c < 0. This makes, 
however, its disagreement with numerical simulations in the region < c < 1 all the 
more surprising. What is it in derivation of Eq. (|l^) that breaks down for c > 0? Or 
are the simulations dominated by finite-size errors and simulations of larger systems 
will eventually agree with Eq. ([l9|)? 

The result for the Ising model is even more interesting. As the theoretical 
predictions are obtained for an Ising model coupled dynamically to the disorder, 
this supports the conjecture put forward in the Introduction about the equivalence 
between annealed and quenched averages over disorder. That is, the only thing 
relevant for the critical behavior of the Ising model are the average fractal properties 
of graphs the spins "see". How the statistical average over graphs is performed, 
quenched or annealed, is not relevant. 

It is also worth noting that we can continuously change the average fractal 
properties of the graphs by changing the embedding dimension D in Eq. (^). This 
allows a continuous interpolation between a flat surface and surfaces corresponding 
to pure gravity. If the prediction of Liouville theory, the KPZ formula, holds for 
all those models, this implies that the critical behavior of the Ising model should 
change continuously in the process. In the language of the renormalization group 
this implies a continuous line of fixed points, rather than isolated points. There 
are well known examples of this; the low-temperature phase of the two-dimensional 
Xy-model or the critical line of the Ashkin- Teller model. But is this statement 
also true for very weak disorder? If we change the fractal dimension infinitesimally, 
from 2 to 2 -|- e, is that enough to change the critical behavior of the Ising model? 
Or, alternatively, does there exist some central charge c' < —5 were the geometrical 
disorder is not strong enough and we always get the Onsager exponents? This point 
deserves further study. 
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One could also look at the examples of weak disorder that have been studied 
recently, for example the site or bond-diluted Ising model, and ask if that kind of 
disorder can also be classified according to some average fractal properties of the 
lattices. And, moreover, if one could observe some kind of universality in the critical 
behavior, depending on the fractal structure, akin to what we have presented in this 
paper. 

In view of how dramatically the critical behavior of the Ising model changes 
on surfaces with such strong disorder, one might ask if such change could be ob- 
served in real physical systems. Possible candidates for such systems could be, for 
example, electrons trapped on the interfaces between two liquids, or on the surface 
of some porous material, were the surfaces had some well defined non-trivial fractal 
structure. As our results indicate, it is only the average geometry of the surfaces 
that is important for the Ising model, not its fluid nature or curvature fluctuations. 
Thus the relative time-scale between the interactions of the particles and the change 
in the geometry should be irrelevant. 

Acknowledgments: The work of G.T. was supported by the Humboldt Founda- 
tion. The work of P.B. was partially supported by KBN grants 2P03 B19609 and 
2P03 B04412. 

References 

[1] A.B. Harris, J. Phys. C7 (1974) 1671. 

[2] Vik.S. Dotsenko and Vl.S. Dotsenko, J. Phys. C15 (1982) 495; Adv. Phys. 32 
(1983) 129; 

V.B. Andreichenko, Vl.S. Dotsenko, W. Selke and J.-S. Wang, Nucl. Phys. 
B344 (1990) 531; Europhys. Lett. 11 (1990) 301; Physics A164 (1990) 221; 
A.L. Talapov and L.N. Shchure, ( |hep-lat/9404002D 
J.K. Kim and A. Patrascioiu, Phys. Rev. Lett. 72 (1994) 2785. 

[3] A.W.W. Ludwig, Nucl. Phys. B330 (1990) 639; Nucl. Phys. B285 (1987) 97. 

[4] V.G. Knizhnik, A.M. Polyakov and A.B. Zamolodchikov, Mod. Phys. Lett. 
A3 (1988) 819. 

[5] C.F. Baillie, K.A. Hawick and D.A. Johnston, Phys. Lett. B328 (1994) 284; 
D.A. Johnston, Phys. Lett. B277 (1992) 405. 

[6] C.F. Baillie, W. Janke and D.A. Johnston, Phys. Lett. B388 (1996) 14; 
Nucl. Phys. 53 (Proc. Suppl.) (1997) 732. 

[7] S. Catterall, G. Thorleifsson, M. Bowick and V. John, Phys. Lett. B354 
(1995) 56; 

[8] J. Ambj0rn, J. Jurkiewicz and Y. Watabiki, Nucl. Phys. B454 (1995) 313. 



21 



[9] N. Ishibashi and H. Kawai, Phys. Lett. B314 (1993) 190; Phys. Lett. B322 
(1994) 67; 

M. Fukuma. N. Ishibashi, H. Kawai and M. Ninomiya, Nucl. Phys. B427 
(1994) 139; 

J. Ambj0m, C.F. Kristjansen and Y. Watabiki, Nucl. Phys. B504 (1997) 555. 

[10] Y. Watabiki, Progress in Theoretical Physics, Suppl. No. 114 (1993) 1; 
H. Kawai, N. Kawamoto, T. Mogami and Y. Watabiki, Phys. Lett. B306 
(1993) 19. 



[11] P. Ginsparg, Matrix Models of 2d- Gravity, ([hep-th/9112013D , Lectures given 
at Trieste Summer School, Trieste, Italy, July 22-25, 1991; 
L. Alvarez-Gaume, Random Surfaces, Statistical Mechanics and String 
Theory, Helv. Phys. Acta 64 (1991) 359. 

[12] F. David, Simplicial Quantum Gravity and Random Lattices, 

( |iep-th/9303127 ), Lectures given at Les Houches Summer School on 
Gravitation and Quantization, Session LVII, Les Houches, France, 1992; 
J. Ambj0rn, Quantization of Geometry, ( |hep-th/941117"9 ). Lectures given at 
Les Houches Summer School on Fluctuating Geometries and Statistical 
Mechanics, Session LXII, Les Houches, France 1994; 

P. Di Francesco, P. Ginsparg and J. Zinn-Justin, Phys. Rep. 254 (1995) 1. 

[13] V.A. Kazakov, Phys. Lett. A119 (1986) 140; 

D. Boulatov and V.A. Kazakov, Phys. Lett. B186 (1987) 379; 
Z. Burda and J. Jurkiewicz, Acta. Phys. Polon. B20 (1989) 949. 

[14] M. Bowick, S.M. Catterah and G. Thorleifsson, Phys. Lett. B391 (1997) 305. 

[15] J. Ambj0rn, G. Thorleifsson and M. Wexler, Nucl. Phys. B439 (1995) 187. 

[16] V.A. Kazakov, I.K. Kostov and A.A. Migdal, Phys. Lett. 157B (1985) 295. 

[17] W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery, Numerical 
Recipes, Cambrigde University Press (1992) 

[18] J. Ambj0rn, Ph. De. Forcrand, F. Koukiou and D. Petritis, Phys. Lett. B197 
(1987) 548. 

[19] S. Jain and S.D. Mathur, Phys. Lett. B286 (1992) 239. 

[20] J. Ambj0rn, S. Jain and G. Thorleifsson, Phys. Lett. B307 (1993) 34; 
J. Ambj0rn and G. Thorleifsson, Phys. Lett. B323 (1994) 7. 

[21] J. Ambj0rn and K.N. Anagnostopoulos, Nucl. Phys. B497 (1997) 445. 

[22] J. Ambj0rn, K.N. Anagnostopoulos, T. Ichihara, L. Jensen, N. Kawamoto, 
Y. Watabiki and K. Yotsuji, Phys. Lett. B497 (1997) 177; Intrinsic Geometric 
Structure of c = —2 Quantum Gravity, Nucl. Phys. B511 (1998) 673. 



22 



[23] P. Bialas, Phys. Lett. B373 (1995) 289. 

P. Bialas, Z. Burda and J. Jurkiewicz, Correlation functions and critical 
behavior on fluctuating geometries, (hep-lat/9710070). 

[24] J. Ambj0rn and Y. Watabiki, Nucl. Phys. B445 (1995) 129. 

[25] J. Ambj0rn, K.N. Anagnostopoulos, U. Magnea and G. Thorleifsson, 
Phys. Lett. B388 (1996) 713. 

[26] M.J. Bowick, V. John and G. Thorleifsson, Phys. Lett. B403 (1997) 197. 

[27] J. Ambj0rn, K.N. Anagnostopoulos, J. Jurkiewicz and C.F. Kristjansen, The 
Concept of Time in 2d Gravity, ([hep-th/980202QD . 

[28] C. Itzykson and J.-M. Drouffe, Statistical Field Theory, Cambridge University 
Press, Cambridge, 1989. 

[29] W. Janke, M. Katoot and R. Villanova, Phys. Lett. B315 (1993) 412; 
Phys. Rev. B49 (1994) 9644; 

W. Janke and R. Villanova, Phys. Lett. A201 (1995) 179. 

[30] V.A. Kazakov, M. Staudacher and T. Wynter, Nucl. Phys. B471 (1996) 309. 

[31] A.M. Ferrenberg and R.H. Swendsen, Phys. Rev. Lett. 63 (1989) 1195; ibid 
61 (1988) 2635. 

[32] K. Binder and A.P. Young, Rev. Mod. Phys. 58 (1986) 801. 



23 



